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Abstract. New test-particle simulations have been per- 
formed to study the secular evolution of the gaseous dis- 
tribution in NGC 4736. We find that the distribution of 
gas clouds can be understood in the frame of perturba- 
tion induced by a nuclear oval potential, in addition to a 
spiral-like potential. Our experiments show that both in- 
ner and outer rings are stable structures located at the in- 
ner Lindblad resonance (ILR) and the outer Lindblad res- 
onance(OLR), respectively, in agreement with the obser- 
vations of NGC 4736. One of our simplified simulations in- 
dicates that both nuclear starburst and orbital resonance 
might be needed for keeping the inner ring stable for a 
long period. We have introduced a symplectic algorithm 
in the orbit integration. Its stronger stability allows us to 
adopt a rather large time step to save computational time. 
Substitution of a viscous force for cloud-cloud interactions 
proves adequate in the case of NGC 4736. 

Key words: galaxies: spiral - galaxies: evolution - galax- 
ies: kinematics and dynamics - galaxies: structure - galax- 
ies: individual (NGC 4736) - methods: numerical 



1. Introduction 

NGC 4736 (M94) is a bright, nearby Sab galaxy 0, which 
is notable for a faint outer Hi ring of radius 4 arcmin 
to 6 arcmin, and an inner bright ring of H II region with 
enhanced gas density at 45 arcsec (Lynds 1974; Beckman 
et al. 1991). Beckman et al. (1991) found that the bulge 
position angle (~ 20°) was very different from that of the 
disc (~ 120°), strong evidence for non-axisymmetry in the 
bulge of NGC 4736, and suggested that triaxiality of the 
bulge played an important dynamical role in fuelling star 
formation. 



Send offprint requests to: Q.-S Gu, e-mail:postcstd@nju. edu.cn 
1 The Hubble distance is still uncertain, we adopt 6.3 Mpc 
in this paper, so that 1 arcmin is equal to 1.8 kpc. 



Discussion of the origin of these rings is separated into 
two groups. On one hand, van der Kruit(1974, 1976) and 
Sanders & Bania(1976) suggested that the ring-like struc- 
tures in NGC 4736 are the observational signature of re- 
cent nuclear explosive events. On the other hand, Schom- 
mer & Sullivan (1976) and Bosma et al. (1977) proposed 
that the rings result from tight winding of spiral arms 
and are tracers for the principle orbital resonance regions. 
Athanassoula et al. ( 1982 ) thought that NGC 4736 is in 
fact the prototype of rings located at Lindblad resonances 
in an oval potential. 

Gerin, Casoli & Combes (1991) have shown that both 
inner H n and outer H I rings could form in a slight barred 
potential. The possible existence of such an oval was first 
discussed by Bosma, van der Hulst & Sullivan (1977) and 
Kormendy (1979) and more recently by Huang, Gu & Su 
(1993). As pointed out by the authors themselves, the 
model of Gerin, Casoli & Combes has two main draw- 
backs. Their outer ring is too wide (cf. their Fig. 10b), 
so that it was in fact associated with both corotation at 
about 4.5 kpc and OLR at about 8.0 kpc. The second 
drawback is that their inner ring disappears after 3 Gyrs. 
The aim of this paper will be to produce a model devoid 
of these drawbacks. 

Here, we present our numerical experiments on gas 
clouds in NGC 4736 based on the test-particle method. 
We replace the process of cloud-cloud collision with the 
viscous force widely used in the study of accretion disks 
and introduce the symplectic algorithm, which has proven 
of strong stability, in our simulations. Both inner and 
outer ring structures turned out to be quite stable after 
4.0 Gyr, and are associated with the inner Lindblad res- 
onance (ILR) and outer Lindblad resonance(OLR) for a 
long duration, respectively. 

In Sect. 2, we describe the construction of the model 
potential of NGC 4736, and the method of replacing the 
process of collision between gas clouds. The computa- 
tional method, namely, the near-symplectic algorithm, is 
described in Sect. 3, and in Sect. 4 we present the simula- 
tion results of the gasdynamics in NGC 4736. Finally, we 
give the discussions and conclusions in Sect. 5 and 6. 
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2. The new simulations 

We assume in our new simulations that the motions of 
the gas clouds are mainly dominated by the following two 
physical processes: the galaxy's gravitational field and the 
mutual collision between gas clouds. Note that we neglect 
the self-gravity of the gas cloud, which may have an im- 
portant effect in some cases. 



2.1. The gravitational field 

Following the basic approach, the galactic potential con- 
sists of two parts, an axisymmetric background and a 
non-axisymmetric perturbing potential. We assume that 
the background potential has four components: a spheri- 
cal bulge, a thick disk, a thin disk, and a spherical halo. 
Both the density distribution of the bulge and of the halo 
take the form 



r > r c . 



(1) 



which implies an underlying bulge gravitational potential 
of the form 



-2.0Gp r 2 (i - £), r < r c 
-Gp rl/6r, r > r c 



(2) 



The halo potential Un has the same form as Ub- The thick 
and thin disks are assumed to be Toomre disks (Toomre, 
1963), and have the form 



Uu(r) 



( r 2 + r 2)l/2' 



(3) 



We take the perturbing potential as the same as that 
used in Roberts & Hausman's simulations (1984). 



U 1 (r,e,t) = U J) (r) 



5 (r 2 + r 2 ) 2 



cos[20-2fi p i + $(r)]. (4) 



where <&(/■), the phase of the maximum perturbation, is 
determined by 



$(r) = 21og(l + (r/r ) j )/(j tan i ), 



(5) 



In the above equations, r c is a constant canonical radius, 
fi p is pattern speed, Md the total mass of the Toomre disk, 
ro a characteristic radius where the perturbing potential 
changes from bar like to spiral-like. The power j determines 
how sharply the shift from barlike to spiral perturbation 
occurs with increasing r, and io is the pitch angle of the 
spiral. In our model, ro, j, and io are set to be l.Okpc, 5, 
and 10°, respectively. 

Now we can construct a mass model of NGC 4736 by 
fitting the observed rotation curve with the model poten- 
tials mentioned above. The rotation curve of NGC 4736 
has the following characteristics: a steep rise from the nu- 
cleus to 0.25 arcmin, a turn-over region from 0.25 to 0.5 



Fig. 1. Rotation curve model of NGC 4736 



Fig. 2. Circular angular velocity ( f2 ) and Lindblad precession 
frequencies ( Q ± n/2 ) vs. radius for the rotation curve model 
of NGC 4736 



arcmin, a constant-velocity region from 0.5 to 4.0 arcmin 
and a decreasing- velocity region from 4.0 to 6.0 arcmin 
(Schommer & Sullivan 1976 ; Bosma et al. 1977). The 
mass parameters, determined this way, of the bulge, thick 
disk, thin disk and halo are 0.005, 0.032, 0.11, and 0.11 
kpc 3 /Myr 2 , respectively, and the constant canonical radii 
for these four components are 0.30, 0.45, 2.50, and 8.20 
kpc, respectively. Fig.l gives the rotation curve model cor- 
responding to the unperturbed gravitational field defined 
by the above parameter values, which is almost in keeping 
with that derived from optical and radio observations. The 
standard Lindblad precession frequencies for this rotation 
curve are shown in Fig. 2. 

The equation of motion of gas clouds in the potential 
U (x, y, t) can be written as 



dx dH dx 

dt dx ' dt 

dy _ dH_ dy 

dt dy ' dt 



dH 
dx ' 
dH 



(6) 
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where, H, the Hamiltonian of the gas cloud, is represented 
by 



H = H(x, y, x, y, t) = -(x 2 + y 2 ) + U(x,y,t), 



(7) 



x and y are the canonical conjugate variables of x and y, 
respectively. 

2.2. The mutual collisions between gas clouds 

The cloud-cloud collision is obviously related to the fol- 
lowing parameters: the gas density, the velocity disper- 
sion, the collision cross-section, the constitution of matter, 
the differential rotation, the included collision angle, etc. 
Therefore, it is difficult to formulate the collision processes 
exactly in numerical simulations. The effect of the mutual 
collision between gas clouds was treated as a simple frac- 
tional reduction of the relative velocities first by Schwarz 
(1981,1984). In our present numerical experiments, we in- 
troduce a dissipative viscous force to simulate the effect of 
collision, expressed by the following formula in the well- 
known a-disk (Shakura and Sunyaev 1973; Pringle 1981), 



aC s T,hr^-, 
dr 



(8) 



where, a is a constant and < a < 1; C s the local sound 
velocity; £ the cloud surface density; f2 the local angular 
velocity; and h the spiral thickness. For the gas clouds in 
a galaxy, C s is just the local velocity dispersion. As a first- 
order approximation, we take v — vq as C s , where v and 
vo is the local velocity and the local unperturbed circular 
rotation speed, the initial number density of gas clouds So 
as £, and the local unperturbed circular rotation f2 as f2. 
So, we have 



F = a(v — v )T, h r 



dfl 
dr 



(9) 



to replace the effect of collision between gas clouds. Ob- 
viously, F is a small quantity of first-order magnitude as 
compared with the unperturbed gravitational force. With 
the addition of F, the motion equation of a single gas 
cloud on the circle of radius r becomes 



dx dH dx dH I n rp 



dx 

dt 
dy 

dt 



dH dx 

dx ' dt 

OH dy 

di) ' dt 



dH 

dx 

— + lirrF 
dv t ^7T' r y , 



(10) 



where F x and F y are the components of F in the directions 
of x and y, respectively. 

3. The computational method 

For the numerical computation of the Hamiltonian system, 
we could adopt traditional methods such as Rungc-Kutta. 
But most of the traditional numerical methods will cause 
a linear variation of the system energy with time, con- 
trary to the theoretical results for Hamiltonian systems. 



It is known that the Hamiltonian flow preserves its sym- 
plcctic structure (Arnold, 1978). Based on this property, 
a symplectic algorithm for Hamiltonian systems has been 
developed (Ruth 1983; Feng 1984). Its obvious advantage 
over the traditional methods is in the aspect of the long- 
term dynamical evolution of Hamiltonian system, i.e. no 
secular change in the system energy. 
For the Hamiltonian system 



dH 

dp ' 

[y dq> 

with Hamiltonian function 
H = H(p,q)=T(p) + V(q), 



(11) 



(12) 



where p,q £ R n are a set of canonical conjugate variables, 
its phase flow is represented by 



(13) 



where the differential operator Dh is defined by the Pois- 
son bracket {•, •} as follows 



(14) 



and exp[(t — to) Da] is the exponential transformation of 
Dh- Let the integration time-step be r, both of the fol- 
lowing step-transition operators 



g\ = exp(TDT)exp(TDy) 



gl — exp(TDy)exp(TDT) 



(15) 



(16) 



are first-order symplectic integrators. It is easy to prove 
that both 



.92 = 9i ° 9i 
and 

92 = 01 ° 01 



(17) 



(18) 



are second-order symplectic algorithms. Higher-order 
symplectic algorithms can be constructed from a certain 
number of second-order symplectic algorithms (Yoshida, 
1990). 

The difference scheme corresponding to ^ 



9k+i — 9k -T 2 g p \p=p k , 

Pk+l = Pk-T^)| 9=9k+i , 

— _1_ r dT{p) | 2 

<Zk+l — <7k+i ' 2 dp lp=Pk+n 



(19) 



where (pk+i, <Zk+i) and (pk, 9k) are the numerical solutions 
of the Hamiltonian flow (13) at time t = (k + 1)t and 
t = kr, respectively, and po = p(to), qo = q(to)- 
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Fig. 3. Time evolution of the gas clouds in our model. The time epochs for a, b, c, d, e, and f are 50 Myr, 200 Myr, 400 Myr, 
1.0 Gyr, 2.0 Gyr and 4.0 Gyr, respectively. 



If we consider the small dissipative effect, then equa- 
tion (11) becomes 

(q= d*L 

\p = - P ^+eF( P , q ), W 

Since e is a small quantity, we could replace the solution 
of the second equation (20) by an Eulcr flow. Hence the 
difference scheme (19) changes to 

(„ , - „, I rdT(p), 
yk+i — * T 2 dp |p=p k , 

< p k+1 = Pk-r[^J) + eF(p,q)}\ p=p ^ q= ^ + ,, (21) 
Qk+i = q k+ i + f^^^pk+i, 

This kind of treatment causes the truncation error of (21) 
to be of order max[0(r 3 ), 0(r 2 e)]. Obviously, we can en- 
sure that the scheme (21) is second-order if e < r, and we 
will use this scheme in our simulation. 

It should be noticed that the "leap-frog" difference 
scheme used in the TV-body simulation is also symplec- 
tic and second-order. But its stable interval is very small 
in comparison with scheme (21). 



4. Results 

At the beginning of each simulation, 10000 gas clouds are 
uniformly distributed in the disk with 10 kpc radius. Each 
cloud is given with a local circular rotation velocity, along 
with a peculiar velocity selected from a two-dimensional 
Gaussian distribution whose one-dimensional dispersion is 
5 km s _1 . 

With the amplitude of the spiral perturbation A = 0.4 
, the pattern speed Sl p = 0.03, and the viscous coefficient 
a = 0.001, the dynamical evolution of gas clouds with 
time is given in Fig. 3. It is clearly shown that when time 
t = 50Myr, Fig. 3a, a faint ring-like structure of radius 
l.lkpc emerges in the central region; with increasing time, 
more and more gas particles move inward and accumulate 
at the ring region. At t = 400Myr, the inner stable ring 
becomes prominent, roughly located at the ILR with a 
radius of 1.1 kpc , see Fig. 3c. At the same time, in the 
outer part of the disk, an outer ring is formed as the result 
of tight winding of spiral arms. It is stable and of low- 
density with radius 7.5 kpc to 10.0 kpc, a location near 
the OLR and beyond. It is interesting to notice that both 
inner and outer ring structures are stable for durations 
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Fig. 4. Time evolution of the gas clouds, to which within the ILR region some radial motions are added when E > E c . The 
time epochs for a, b and c are 500 Myr, 1.0 Gyr, and 3.6 Gyr, respectively. 



of 100 Myr to 800 Myr, fixed onto the ILR and OLR, 
respectively, in agreement with the observations of NGC 
4736. 

After 1.0 Gyr, the inner ring moves inward further, to 
a radius of about 0.7kpc, see Fig. 3d, and stays there till 

4.0 Gyr, see Fig. 3f. Meanwhile, the outer ring appears 
narrower, but remains in the OLR region. 

5. Discussion 

5.1. Validity 

The results obtained in Fig 3. are conditional on the ab- 
sence of other complicated processes during the cloud- 
cloud collision, e.g., the formation of giant molecular 
clouds and subsequent instabilities, star formation or mas- 
sive star bursts, and such follow-up processes as heating 
by stars, stellar wind, protostellar jets and poststellar ex- 
plosions. Indeed, such processes could be triggered when 
the surface density of gas clouds becomes high enough 
(Elmegreen 1994), and the motions and thermodynamic 
state of the gas clouds could be affected (Krugel & Tu- 
tukov 1993). Obviously, our present simulation code will 
not be valid in such case. 

We give a very simplified example. Bcckman et al ( 
1991 ) suggest that starburst activities around the nu- 
cleus of NGC 4736 can produce an outflow motion of 
~ 30km s _1 , as observed by van der Kruit ( 1976 ), Thus 
we add radial motions (~ 30km s^ 1 ) to some of the gas 
clouds inside the ILR when the surface density £ is higher 
than a critical value E c = a 3 ^ G (Toomre 1964, Cowie 
1981, Kennicutt 1989). This might be a way to understand 
the possible influence of star formation on the distribution 
of gas clouds in the ILR region. The simulation obtained 
in this case is shown in Fig 4. It is interesting to find that 
the inner ring remains stable at the ILR location of about 

1.1 — 1.2 kpc even after 3.6 Gyr. Although this approach 



to simulating the effect of star formation is very crude, it 
does give a clue that the debate on the origin of the ring 
structures in NGC 4736 mentioned in Sect. 1 could be 
reconciled in that both nuclear explosion, an event proba- 
bly similar to massive starburst, and orbital resonance are 
needed for keeping the inner ring stable for a long period 
of time. More work on this matter is under consideration. 

5.2. Modelling of the cloud-cloud interactions 

In comparison with others, and beyond the potential 
adopted for NGC 4736, we have modified both the orbit in- 
tegrator and the modelling of the cloud-cloud interactions. 
To see whether the latter is the basic factor improving our 
results, as suggested by Lia Athanassoula ( 1995, private 
communication ), we have repeated one of Schwarz's sim- 
ulations ( 1984 ), using exactly the same potential as he, 
but our model for cloud-cloud interactions. The results are 
shown in Fig 5 with the viscous coefficient a = 0.0001. It 
is interesting to note that the depopulation of the L4, L5 
Lagrangian points is much more rapid than Schwarz's re- 
sults, shown in his Fig 6 ( 1984 ). The outer ring forms 
after four bar rotations, and becomes almost circular and 
stable after six bar rotations. The same structures formed 
in Schwarz's simulation after twenty bar rotations. It ap- 
pears that the basic difference between Schwarz's simula- 
tion and ours is the time evolution and the time scale for 
stabilization. Thus, we believe that the different modelling 
of the cloud-cloud interactions might be the dominant fac- 
tor responsible for the improvement in our results. 

5.3. The approximation 

An important approximation in our numerical experi- 
ments was to keep the surface density, £, in eqn.(8) con- 
stant during the time evolution of gas clouds. Apparently, 
it is far from being the CclSG, clS the simulation indicated. 
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Fig. 5. Time evolution of the gas clouds, using Schwarz's potential but our model for cloud-cloud interactions. The time epochs 
in bar rotations are 2, 4, 6, 8, 10, and 20 for frames a, b, c, d, e, and f, respectively. 



Fig. 6. Time evolution of the gas clouds, same parameters as those in Fig. 3 were used except that the variation of E was 
tracked at each time step. The time epochs for a, b and c are 50 Myr, 400 Myr and 1.0 Gyr, respectively. 
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We have examined the effect of this approximation. 
Instead of keeping constant surface density, So, a more 
time-consuming simulation has been performed, in which 
the variation of surface density was tracked for each step. 
One such simulation is shown in Fig 6 . As compared with 
those in Fig. 3, one can find major differences in the inner 
ring, which moves closer to the nucleus and gathers many 
more gas clouds. On the other hand, the approximation 
has little effect on the outer ring structure. 

Though tracking the variation of £ is more realistic, 
that more gas clouds assemble near the galactic nucleus 
may not be true. As pointed out in Sect. 5.1, when the gas 
density exceeds the critical value E c , the process of star 
formation or even starburst must be considered instead. 
The resulting superwind might prohibit the gas clouds 
from moving inward, as our simple simulation in Fig 4 
indicated. So, it might not be reasonable to track the vari- 
ation of £ without considering the processes of intense star 
formation and the subsequent superwind. In other words, 
the approximation made in the present simulation might 
not be too bad. 

6. Conclusions 

In this paper, we present a new "sticky particle" model 
of NGC 4736, which has three improvements over previ- 
ous models for this galaxy: (1) The orbit integrator; (2) 
the potential adopted for the galaxy; and (3) the different 
modelling of cloud-cloud interactions. The main advan- 
tages are as follows: (1) The numerical experiments show 
that both inner and outer rings in NGC 4736 are sta- 
ble structures located at the ILR and OLR, respectively, 
more in agreement with the observations. The simulations 
overcome the drawbacks described in Sect. 1; (2) The sym- 
plectic algorithm has an advantage over the usual leapfrog 
method in the time step taken, which could save a lot of 
computational time during the simulations; (3) The major 
improvement is due to the different modelling of the cloud- 
cloud interactions with the viscous force as suggested in 
Sect. 5.2. 

Our very simplified simulation gives us a clue that both 
nuclear explosions caused by starburst activities and orbitl 
resonance might be needed to keep the inner ring stable. A 
more realistic way to introduce the process of star forma- 
tion as well as the cloud-cloud collision into our simulation 
is under consideration. 

Acknowledgements. We would like to thank the referee, Dr. E. 
Athanassoula, for her careful reading our manuscript and her 
instructive suggestions, that significantly improve our analy- 
ses in this paper. We would also like to thank the anonymous 
language corrector for improving and correcting our English 
presentation. This research was supported by grants from Na- 
tional Science and Technology Commission and National Nat- 
ural Science Foundation of China. 



Arnold V.I., 1978, Mathematical Method of Classical Mechan- 
ics, Springer, New- York. 

Athanassoula, E., Bosma, A., Creze, M., & Schwarz, M.P., 
1982, A&A, 107, 101 

Beckman J. E., Varela A. M., Munoz-Tunon C, Vilchez J. M., 
Cepa J., 1991 A&A, 245, 436 

Bosma A., van der Hulst J. M., Sullivan III W. T., 1977, A&A, 
57, 373 

Cowie L. L., 1981, ApJ, 245, 66 
Elmegreen B. G., 1994, ApJ 425, L73 

Feng K. 1985 in Proc. 1984 Beijing Symp. Diff. Geom. and 

Diff. Eq. Ed. K. Feng, Science Press, Beijing, 42 
Gerin M., Casoli F., Combes F., 1991, A&A, 251, 32 
Lynds B. T., 1974, ApJSS, 28, 267 

Huang J. H., Gu Q. S., Su H. J., 1993, in Mass-induced Ac- 
tivities in Galaxies, ed I. Shlosman, Cambridge University 
Press, P. 119 

Kennicutt R. C, 1989, ApJ, 344, 685 

Kormendy J., 1979, ApJ, 227, 714 

Kmgcl E., Tutukov A. V., 1993, A&A, 275, 416 

Pringle J. E., 1981, ARA&A, 19,137 

Roberts W. W., Huntley J. M., van Albada G.D., 1979, ApJ, 
233, 67 

Roberts W. W., Hausman M. A., 1984, ApJ, 277, 744 
Ruth R. D., 1983, IEEE, Trans, on Nuc. Sci., NS-30, 2669 
Sanders R. H., Bania T. M., 1976, ApJ, 204, 341 
Schommer R. A., Sullivan III W. T., 1976, Astrophys. Lett., 
17, 191 

Schwarz M. P., 1981, ApJ, 247, 77 

Schwarz M. P., 1984, MNRAS, 209, 93 

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337 

Toomre A., 1963, ApJ, 138, 385 

Toomre A., 1964, ApJ, 139, 1217 

van der Kruit P. C, 1974, ApJ, 188, 3 

van der Kruit P. C, 1976, A&A, 52, 85 

Yoshida H., 1990, Phys. Letters A, 150, 262 



References 



This article was processed by the author using Springer- Verlag 
WT E X A&A style file L-AA version 3. 



